home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 1998 November / IRIX 6.5.2 Base Documentation November 1998.img / usr / share / catman / p_man / cat3 / complib / dlatps.z / dlatps
Text File  |  1998-10-30  |  7KB  |  199 lines
  1.  
  2.  
  3.  
  4. DDDDLLLLAAAATTTTPPPPSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DLATPS - solve one of the triangular systems   A *x = s*b or A'*x = s*b
  10.      with scaling to prevent overflow, where A is an upper or lower triangular
  11.      matrix stored in packed form
  12.  
  13. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  14.      SUBROUTINE DLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM,
  15.                         INFO )
  16.  
  17.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  18.  
  19.          INTEGER        INFO, N
  20.  
  21.          DOUBLE         PRECISION SCALE
  22.  
  23.          DOUBLE         PRECISION AP( * ), CNORM( * ), X( * )
  24.  
  25. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  26.      DLATPS solves one of the triangular systems transpose of A, x and b are
  27.      n-element vectors, and s is a scaling factor, usually less than or equal
  28.      to 1, chosen so that the components of x will be less than the overflow
  29.      threshold.  If the unscaled problem will not cause overflow, the Level 2
  30.      BLAS routine DTPSV is called. If the matrix A is singular (A(j,j) = 0 for
  31.      some j), then s is set to 0 and a non-trivial solution to A*x = 0 is
  32.      returned.
  33.  
  34.  
  35. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  36.      UPLO    (input) CHARACTER*1
  37.              Specifies whether the matrix A is upper or lower triangular.  =
  38.              'U':  Upper triangular
  39.              = 'L':  Lower triangular
  40.  
  41.      TRANS   (input) CHARACTER*1
  42.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  43.              (No transpose)
  44.              = 'T':  Solve A'* x = s*b  (Transpose)
  45.              = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
  46.  
  47.      DIAG    (input) CHARACTER*1
  48.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  49.              Non-unit triangular
  50.              = 'U':  Unit triangular
  51.  
  52.      NORMIN  (input) CHARACTER*1
  53.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  54.              contains the column norms on entry
  55.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  56.              computed and stored in CNORM.
  57.  
  58.  
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDLLLLAAAATTTTPPPPSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      N       (input) INTEGER
  75.              The order of the matrix A.  N >= 0.
  76.  
  77.      AP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
  78.              The upper or lower triangular matrix A, packed columnwise in a
  79.              linear array.  The j-th column of A is stored in the array AP as
  80.              follows:  if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  81.              if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  82.  
  83.      X       (input/output) DOUBLE PRECISION array, dimension (N)
  84.              On entry, the right hand side b of the triangular system.  On
  85.              exit, X is overwritten by the solution vector x.
  86.  
  87.      SCALE   (output) DOUBLE PRECISION
  88.              The scaling factor s for the triangular system A * x = s*b  or
  89.              A'* x = s*b.  If SCALE = 0, the matrix A is singular or badly
  90.              scaled, and the vector x is an exact or approximate solution to
  91.              A*x = 0.
  92.  
  93.      CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
  94.  
  95.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  96.              the norm of the off-diagonal part of the j-th column of A.  If
  97.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  98.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  99.              greater than or equal to the 1-norm.
  100.  
  101.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  102.              the 1-norm of the offdiagonal part of the j-th column of A.
  103.  
  104.      INFO    (output) INTEGER
  105.              = 0:  successful exit
  106.              < 0:  if INFO = -k, the k-th argument had an illegal value
  107.  
  108. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  109.      A rough bound on x is computed; if that is less than overflow, DTPSV is
  110.      called, otherwise, specific code is used which checks for possible
  111.      overflow or divide-by-zero at every operation.
  112.  
  113.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  114.      A is lower triangular is
  115.  
  116.           x[1:n] := b[1:n]
  117.           for j = 1, ..., n
  118.                x(j) := x(j) / A(j,j)
  119.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  120.           end
  121.  
  122.      Define bounds on the components of x after j iterations of the loop:
  123.         M(j) = bound on x[1:j]
  124.         G(j) = bound on x[j+1:n]
  125.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDLLLLAAAATTTTPPPPSSSS((((3333FFFF))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      Then for iteration j+1 we have
  141.         M(j+1) <= G(j) / | A(j+1,j+1) |
  142.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  143.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  144.  
  145.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  146.      j+1 of A, not counting the diagonal.  Hence
  147.  
  148.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  149.                      1<=i<=j
  150.      and
  151.  
  152.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  153.                                       1<=i< j
  154.  
  155.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTPSV if the
  156.      reciprocal of the largest M(j), j=1,..,n, is larger than
  157.      max(underflow, 1/overflow).
  158.  
  159.      The bound on x(j) is also used to determine when a step in the columnwise
  160.      method can be performed without fear of overflow.  If the computed bound
  161.      is greater than a large constant, x is scaled to prevent overflow, but if
  162.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  163.      non-trivial solution to A*x = 0 is found.
  164.  
  165.      Similarly, a row-wise scheme is used to solve A'*x = b.  The basic
  166.      algorithm for A upper triangular is
  167.  
  168.           for j = 1, ..., n
  169.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  170.           end
  171.  
  172.      We simultaneously compute two bounds
  173.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  174.           M(j) = bound on x(i), 1<=i<=j
  175.  
  176.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  177.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  178.      bound on x(j) is
  179.  
  180.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  181.  
  182.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  183.                          1<=i<=j
  184.  
  185.      and we can safely call DTPSV if 1/M(n) and 1/G(n) are both greater than
  186.      max(underflow, 1/overflow).
  187.  
  188.  
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.